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Abstract 

High-efficiency  cooling  systems  are  key  points  in  PEMFC  transport  applications,  as  the  volume  constraints  force  the  reduction  of  the  stack 
size  while  increasing  the  power  density.  Moreover,  to  ensure  an  optimal  electrochemical  reaction  over  the  whole  polymer  membrane  surface  and 
hence  a  maximum  efficiency,  the  temperature  field  in  the  cell  must  be  uniform  and  stay  in  a  narrow  range,  around  80-90  °C.  This  study  focuses 
on  improving  the  thermal  performance  of  heat-exchangers  integrated  in  the  bipolar  plates  of  PEMFCs.  The  current  design  of  the  heat-exchangers 
in  these  applications  is  quite  simple;  cooling  liquid  (water)  flows  in  straight  channels  or  serpentines  in  the  rear  of  the  plates.  The  flow  regime  is 
laminar  with  a  Reynolds  number  around  200.  In  order  to  enhance  convective  heat  transfer,  we  propose  here  to  promote  three-dimensional  flow 
inside  cooling  channels  using  a  novel  channel  geometry  that  generates  chaotic  advection  flow.  However,  to  limit  the  size  and  the  electric  resistance 
of  the  bipolar  plates,  the  thickness  must  be  severely  limited.  This  work  concentrates  on  developing  and  characterizing  heat-exchangers  that  can  be 
easily  reduced  in  size  while  preserving  high  thermal  performance. 

©  2005  Elsevier  B.V.  All  rights  reserved. 
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1.  Introduction 

A  PEMFC  system  produces  as  much  electrical  power  as  ther¬ 
mal  power.  In  transport  applications  in  which  the  power  density 
of  the  fuel  cell  is  relatively  high,  much  thermal  energy  must  be 
evacuated  from  the  system.  For  a  fuel  cell  stack,  current  heat- 
exchanger  technology  involves  carving  channels  of  a  depth  equal 
to  half  (or  one-third)  of  the  bipolar  plate  in  which  water  flows, 
as  shown  in  Fig.  1.  The  conventional  heat-exchanger  design  in 
the  bipolar  plate  is  a  simple  network  of  parallel  straight  chan¬ 
nels.  The  objective  for  future  generations  of  PEMFC  stacks, 
however,  is  to  increase  the  working  temperature  from  80-90  to 
120  °C,  thus,  the  heat  generated  will  increase  greatly,  current 
heat-exchanger  design  will  rapidly  reach  its  limits  and  a  new 
design  must  be  found.  This  is  the  objective  of  the  present  paper. 

Among  the  parameters  that  influence  the  PEMFC  perfor¬ 
mance,  the  temperature  homogeneity  over  the  membrane  surface 
is  crucial.  The  temperature  must  be  high  enough  (typically  over 
80  °C)  to  ensure  a  good  electrochemical  reaction  but  low  enough 
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to  prevent  destroying  the  polymer  membrane.  Another  important 
parameter  is  the  electric  resistance  of  the  collector,  it  must  be  as 
low  as  possible  in  order  to  minimize  ohmic  losses,  so  that  the 
heat-exchanger  located  inside  the  collector  must  be  very  thin. 

Hydrodynamic  conditions  in  the  heat-exchanger  in  bipolar 
plates  are  such  that  the  Reynolds  number  of  the  flow  is  around 
200,  which  corresponds  to  the  laminar  regime.  The  anticipated 
future  size  reduction  of  PEMFC  will  reduce  this  Reynolds  num¬ 
ber.  It  is  well  known  that  the  laminar  regime  is  particularly 
unfavorable  to  high  convective  heat  transfer.  With  these  tech¬ 
nological  constraints  in  mind  (temperature  homogeneity,  low 
thickness  and  laminar  regime),  we  propose  in  this  paper  a  type 
of  geometry  that  can  greatly  improve  current  heat- exchanger 
performance  for  PEMFC  cooling  applications. 

Section  1  of  this  paper  presents  the  different  geometries  con¬ 
sidered  and  Section  2  reports  the  numerical  procedure  used  to 
evaluate  the  thermal  performance  of  the  different  geometries. 
Section  3  presents  the  results  and  discussions  and  Section  4  is 
devoted  to  the  conclusions  and  future  work. 

2.  Studied  geometries 

Instead  of  considering  the  whole  heat-exchanger  with  all  its 
cooling  channels,  we  evaluate  the  thermal  performance  of  a  sin¬ 
gle  channel  in  the  network.  The  next  step  (yet  to  be  done)  will 
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Fig.  1.  Two  elementary  cells  of  a  PEMFC  stack  with  a  plate  heat-exchanger  in 
the  bipolar  plate. 


be  to  resolve  the  heat  transfer  problem  for  the  entire  network 
of  several  parallel  channels  separated  by  conductive  materi¬ 
als  in  order  to  evaluate  the  thermal  performance  of  the  whole 
heat-exchanger.  A  study  must  be  undertaken  to  find  the  best 
composition  among  the  channels,  this  point  is  not  addressed 
here. 

Three  different  single  channel  geometries  were  considered, 
see  Figs.  2-4.  Fig.  2  shows  the  classical  case  of  the  straight 
channel.  The  cross-section  is  rectangular  with  an  aspect  ratio 
of  2,  the  channel  height  is  1  mm,  the  hydraulic  diameter  D h  is 
1.33  mm.  These  values  are  typical  of  actual  PEMFC  stacks.  The 
channel  length  is  8  cm. 

Fig.  3  shows  a  zigzag  channel  with  identical  cross-section  and 
hydraulic  diameter.  This  channel  is  a  succession  of  alternating 


Fig.  3.  Zigzag  channel. 


90°  bends,  so  that  the  geometry  is  periodic  in  space  with  a  period 
of  two  alternating  bends.  The  unfold  length  of  a  period  is  equal  to 
1 8  mm.  The  ratio  of  the  curvature  radius  to  the  hydraulic  diame¬ 
ter  in  each  bend  is  5.5.  With  the  considered  Reynolds  number,  Re 


about  200,  the  resulting  Dean  number  is:  D e  =  =  96. 

Under  this  hydrodynamic  condition,  a  secondary  flow  consisting 
of  two  counter-rotating  vortices,  called  Dean  vortices,  is  gener¬ 
ated  in  each  bend  under  the  action  of  centrifugal  forces.  This  flow 
topology  is  known  to  lead  to  increased  convective  heat  transfer 
over  a  straight  channel  [1,2]. 

Fig.  4  shows  the  so-called  “C-shaped”  geometry,  first  intro¬ 
duced  by  Liu  et  al.  [3].  This  is  a  three-dimensional  geometry 
contrary  to  the  two  previous  ones.  The  term  “C-shaped”  is  used 
because  of  the  channel  shape  in  the  upper  plane.  Liquid  flows  in 
a  straight  channel  in  the  lower  plane.  The  channel  cross-section 
here  is  rectangular  with  =  1 .33  mm,  not  square  as  in  Ref.  [3]. 
The  basic  geometry  constituted  of  a  “C-shaped”  section  in  the 
upper  plane  and  a  straight  section  in  the  lower  plane  is  period¬ 
ically  repeated  in  the  longitudinal  direction.  The  unfold  length 
over  one  period  is  equal  to  the  one  for  the  zigzag  geometry,  i.e. 
18  mm. 

Using  a  Poincare  map,  Ref.  [4]  shows  that  fluid  particles  flow¬ 
ing  in  this  geometry  may  have  chaotic  trajectories.  Some  other 
parts  of  the  flow  do  not  exhibit  this  behavior  because  the  flow  is 
regular.  The  presence  of  chaotic  trajectories  in  the  flow  greatly 
enhances  mixing  [5].  It  is  believed  that  this  special  feature  will 
induce  a  large  heat  transfer  coefficient,  but  this  has  not  been 
proved  so  far. 


Fig.  2.  Straight  channel. 


Fig.  4.  “C-shaped”  channel. 
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3.  Numerical  procedure 

In  order  to  evaluate  the  thermal  performance  of  each  geom¬ 
etry,  numerical  computations  were  performed  using  the  com¬ 
mercial  CFD  code  Fluent®.  A  uniform  surface  heat  flux  is 
imposed  on  a  single  channel  surface.  Current  PEMFC  stacks 
generate  about  5000  Wm-2;  in  anticipation  of  future  power 
density  increases,  we  used  cp=  10,000 Wm-2  in  our  calcula¬ 
tions.  The  inlet  section  of  the  channel  is  hydrodynamically  fully 
developed.  An  approached  fully  developed  velocity  profile  for 
rectangular  ducts  is  imposed  [6] .  The  forced  convection  regime 
is  imposed  and  the  Navier- Stokes  equations  are  solved  under 
laminar  assumptions. 

A  parametric  study  was  performed  to  identify  the  grid  mesh 
that  assures  the  capture  of  the  thermal  and  velocity  gradients 
near  the  walls.  A  spatial  resolution  of  40  x  40  meshes  in  the 
cross-section  and  40  meshes  in  2  mm  in  the  ^-direction  was  found 
fully  satisfactory.  The  convergence  of  computations  was  stopped 
when  the  residues  were  less  than  10-7. 

4.  Results  and  discussion 

The  thermal  performance  of  the  single  channel  was  first 
characterized  in  terms  of  the  evolution  along  the  curvilinear 
coordinate  of  the  non-dimensional  convective  heat  transfer  coef¬ 
ficient,  Nu ,  defined  as: 

Nu  =  - - - — 

Tm  —  Tw  X 

where  Tm  is  the  bulk  mean  fluid  temperature  over  the  cross- 
sectional  area  of  the  channel,  Tw  the  channel  wall  temperature 
and  A  is  the  thermal  conductivity  of  the  flowing  fluid  (here  water, 
A  =  0.6Wm-1  K-1).  Fig.  5  presents  the  evolution  of  the  Nus- 
selt  number  with  the  curvilinear  coordinate  for  these  channel 
geometries.  The  Reynolds  number  of  the  flow  is  200. 

For  the  straight  channel,  a  classical  power-law  decrease  is 
found  with  an  asymptotic  value  of  Nu  =  3.00.  This  value  is  con¬ 
sistent  with  the  results  in  Ref.  [7]  and  can  be  taken  as  a  validation 
of  the  CFD  computations. 

A  significant  decrease  in  the  Nusselt  number  is  observed  in 
the  entrance  region  of  the  zigzag  channel  corresponding  to  a 


Fig.  5.  Evolution  of  the  Nusselt  number  with  the  curvilinear  coordinate  for  the 
three  channel  geometries. 


Fig.  6.  Secondary  flow  at  the  exit  of  the  zigzag  channel. 


distance  of  about  0.005  m  beyond  the  channel  inlet.  After  the 
first  period,  the  variation  of  the  Nusselt  number  is  periodic  with 
a  mean  value  of  8.2,  a  period  equivalent  to  two  alternating  bends, 
and  a  maximum  relative  variation  of  about  1.1.  These  results  are 
in  good  agreement  with  those  in  Ref.  [1].  Heat  transfer  is  higher 
in  a  curved  tube  than  in  a  straight  tube  because  of  the  presence 
of  Dean  vortices.  Fig.  6  shows  the  secondary  flow  at  the  exit  of 
the  zigzag  geometry. 

Dean  vortices  transport  cold  particles  from  the  center  of  the 
duct  to  the  hot  regions  close  to  the  wall.  The  secondary  flow 
generates  Nusselt  number  oscillations,  though  the  amplitude 
and  frequency  of  these  oscillations  depend  on  many  param¬ 
eters.  Unlike  a  single  curved  duct,  where  the  oscillations  are 
attenuated  in  the  asymptotic  zone,  in  the  zigzag  geometry,  the 
oscillations  are  kept  constant  by  the  zigzag  perturbation.  This 
geometry  presents  a  clear  improvement  in  thermal  performance 
over  the  straight  channel,  the  mean  Nusselt  number  is  about  2.5 
times  higher. 

After  the  first  period  in  the  C-shaped  channel,  a  periodic  vari¬ 
ation  of  the  Nusselt  number  is  observed  around  a  mean  value 
of  20.4.  The  period  is  equal  to  the  spatial  period  of  the  channel 
and  the  maximum  relative  variation  of  the  Nusselt  number  is 
about  10.  The  mean  value  remains  constant  over  the  computa¬ 
tional  domain  and  is  much  higher  than  that  for  the  two  other 
configurations.  Even  though  the  zigzag  configuration  increases 
heat  transfer,  Dean  vortices  divide  the  cross-section  into  two 
zones  where  the  isotherms  form  closed  curves.  Fluid  particles 
in  the  core  of  Dean  vortices  are  prevented  from  approaching 
the  hot  walls.  The  flow  is  regular  and  no  chaotic  region  exists. 
However,  in  the  C-shaped  geometry,  the  existence  of  chaotic 
advection  regime  can  be  shown.  Actually,  the  chaotic  advection 
is  present  in  a  flow  when  at  least  one  of  the  following  features 
exists  [5]: 

-  strong  stretching  and  folding  of  material  lines  and  surfaces; 

-  sensitivity  to  initial  conditions; 

-  production  of  transverse  homoclinic  and  heteroclinic  intersec¬ 
tion. 

The  second  feature  was  emphasized  in  the  C-shaped  geome¬ 
try.  Figs.  7  and  8  show  the  evolution  of  two-particle  trajectories, 
initially  close  to  each  other  and  injected  at  the  center  of  the 
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inlet  area  of  the  zigzag  and  C- shaped  channels.  The  distance 
separating  the  two  particles  is  equal  to  the  mesh  size.  In  the 
zigzag  channel,  the  flow  is  regular  and  the  particle  trajectories 
do  not  diverge.  In  the  C-shaped  channel,  the  particle  trajecto¬ 
ries  can  diverge  rapidly.  The  distance  separating  two-particle 
trajectories  was  computed  in  several  cross-section  of  the  chan¬ 
nel.  The  distance  increases  exponentially,  a  clear  signature  of 
chaotic  advection,  as  in  Ref.  [4].  Consequently,  in  the  C-shaped 
channel,  the  cold  fluid  in  the  center  is  no  longer  contained  in  the 
cold  region.  More  cold  particles  of  fluid  visit  hot  regions  close 
to  the  wall  and  global  heat  transfer  is  improved.  Heat  transfer 
is  increased  due  to  the  kinematic  effect  and  the  mean  Nusselt 
number  retains  the  same  value  all  along  the  geometry.  These 
results  are  in  good  agreement  with  those  in  Ref.  [8]  for  a  similar 
type  of  chaotic  geometry.  Figs.  9  and  10  present  the  tempera- 


Fig.  9.  Temperature  profiles  at  0.05  m  from  the  inlet  for  the  three  geometries 
(horizontal  direction). 
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Fig.  10.  Temperature  profiles  at  0.05  m  from  the  inlet  for  the  three  geometries 
(vertical  direction). 

ture  profile  for  the  three  geometries  at  0.05  m  from  the  inlet.  In 
both  horizontal  and  vertical  directions,  the  temperature  profile 
is  more  uniform  for  the  C-shaped  channel.  The  flatter  tempera¬ 
ture  distribution  in  the  C-shaped  channel  is  attributed  to  chaotic 
trajectories  of  fluid  particles  due  to  the  generation  of  chaotic 
advection  regime. 

Two  points  about  the  C-shaped  geometry  remain  to  be 
addressed:  first,  the  pressure  loss,  which  is  undoubtedly  higher 
than  in  the  two  other  geometries,  and  secondly,  the  vertical  exten¬ 
sion,  which  doubles  the  space  needed  for  the  C-shaped  geometry 
in  the  bipolar  plate.  The  pressure  loss  is  computed  for  the  three 
geometries  for  a  Reynolds  number  of  200  and  is  expressed  in 
term  of  the  friction  factor  defined  as: 

d P\  Dh 

d  s)  1/2  pUl 

where  p  is  the  fluid  density  and  dp/ds  is  the  local  pressure  gra¬ 
dient  along  the  curvilinear  coordinate  of  the  channel.  As  this 
parameter  is  Reynolds  dependent,  it  is  often  preferable  to  cal¬ 
culate  the  product /x  Re.  Fig.  11  presents  the  evolution  of  the 
term  fRe  with  the  curvilinear  coordinate  for  the  three  geometries 
for  a  Reynolds  number  of  200.  For  the  straight  channel,  as  an 
approached  fully  developed  velocity  profile  has  been  imposed 
at  the  inlet  boundary  of  the  computational  domain,  the  quantity 
fRe  has  an  asymptotic  behavior  toward  a  value  of  62.  This  value 


Fig.  11.  Evolution  of  the  term  fRe  with  the  curvilinear  coordinate  for  the  three 
geometries. 
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Fig.  12.  Evolution  of  the  ratio  of  the  friction  factor  to  the  Nusselt  number  with 
the  curvilinear  coordinate  for  the  three  geometries. 

is  to  be  compared  with  the  theoretical  value  of  fRe  for  a  fully 
developed  pipe  flow  which  is  62.2.  This  result  gives  another 
validation  for  the  present  computations.  The  presence  of  bends 
in  the  zigzag  geometry  increases  the  pressure  loss  over  that  in 
a  straight  geometry.  After  the  establishment  region,  the  mean 
value  of  the  quantity  fRe  for  the  zigzag  geometry  is  about  80. 
The  sharp  angles  in  the  C-shaped  geometry  lead  to  a  mean  value 
of  about  160.  This  value  is  2.6  times  higher  than  that  in  the 
straight  geometry. 

The  question  arises  whether  the  higher  Nusselt  number  for 
the  C-shaped  geometry  is  only  due  to  the  higher  friction  coeffi¬ 
cient.  In  other  words,  is  the  heat  transfer  more  intensified  than 
the  pressure  loss  in  the  C-shaped  geometry?  To  answer  this  ques¬ 
tion,  the  ratio  of  the  friction  factor  to  the  Nusselt  number,  , 
is  introduced.  The  geometry  that  produced  the  most  efficient 
convective  heat  transfer  is  the  one  for  which  the  maximum  Nus¬ 
selt  number  is  coupled  with  the  minimum  friction  factor,  i.e. 
the  smallest  ratio  Fig.  12  shows  the  evolution  of  this  ratio 
along  the  three  channels.  In  the  first  period  of  the  zigzag  and  the 
C-shaped  geometry,  i.e.  for  s  <0.0 18  m,  the  ratio  is  compa¬ 
rable  for  the  three  considered  geometries.  However,  for  larger 
curvilinear  coordinate,  the  C-shaped  geometry  presents  a  sig¬ 
nificant  improvement  in  terms  of  effectiveness  compared  to  the 
two  others;  the  mean  value  of  for  the  C-shaped  geometry  is 
less  than  half  of  that  for  the  straight  geometry.  We  thus  see  that 
with  the  C-shaped  channel,  more  heat  transfer  is  promoted  than 
friction  is  produced. 

It  may  be  argued  that  the  C-shaped  geometry  has  a  signifi¬ 
cant  disadvantage  in  heat-exchangers  for  PEMFC  applications 
because  of  its  height;  in  the  present  study,  it  is  actually  twice  the 
height  of  the  straight  or  zigzag  geometries.  The  first  solution  is  to 
halve  the  size  of  the  C-shaped  geometry  and  keep  the  Reynolds 
number  constant  by  doubling  the  mean  velocity.  If  this  solution 
is  judged  unsatisfactory,  the  Reynolds  number  can  be  reduced 
as  much  as  the  size.  We  can  postulate  that  even  if  the  Reynolds 
number  is  significantly  reduced,  the  Nusselt  number  will  still  be 


higher  than  in  the  current  straight  channel  because  very  high  con¬ 
vective  heat  transfer  has  been  obtained  for  a  Reynolds  number 
of  200  (about  six  times  higher  than  that  for  the  straight  channel). 

5.  Conclusions  and  future  work 

This  numerical  study  considered  two  channel  geometries  as 
alternatives  to  the  one  used  in  present  heat-exchangers  integrated 
in  the  bipolar  plate  of  a  PEMFC  stack.  The  first  alternative,  the 
zigzag  channel,  leads  to  convective  heat  transfer  coefficients 
about  twice  that  of  the  straight  channel,  an  improvement  clas¬ 
sically  attributed  to  the  presence  of  Dean  roll  cells  in  the  flow. 
The  second  alternative  is  the  C-shaped  geometry,  for  which  the 
convective  heat  transfer  coefficient  is  about  six  times  that  of  the 
straight  channel.  This  high  performance  is  explained  as  a  con¬ 
sequence  of  chaotic  regions  in  the  flow.  The  chaotic  behavior 
of  fluid  trajectories  has  been  proven  numerically.  The  pressure 
loss  along  the  C-shaped  channel  is  greater  than  in  the  straight 
channel.  However,  the  ratio  of  the  friction  factor  to  the  Nusselt 
number  is  lower,  showing  that  the  heat  transfer  enhancement  is 
higher  than  the  pressure  loss  increase. 

Future  work  will  focus  on  the  thermal  performance  of  the 
C-shaped  channel  at  low  Reynolds  number  and  on  finding  alter¬ 
native  geometries  to  the  C-shaped  geometry  in  which  large 
chaotic  regions  and  low  pressure  losses  are  promoted. 
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